Preliminary test of different forecasting methods applied to simulated data

Load and tidy up the time series library

library(EWSmethods)
#ewsnet_init(envname = "EWSNET_env", auto = FALSE)
library(earlywarnings)
library(dplyr)
library(tidyverse)
library(patchwork)
library(cowplot, warn.conflicts=FALSE)
library(MuMIn) %>% suppressMessages()
library(chngpt) %>% suppressMessages()
library(ggtext)
library(pheatmap)
library(ggstatsplot)
library(asdetect) # devtools::install_github("caboulton/asdetect")
library(ggradar) # devtools::install_github("ricardo-bion/ggradar")
library(pracma)
library(MatrixModels)
setwd("~/Documents/biodicee/fish@risk")
#setwd("~/Documents/BIODICEE/fishATrisk")
#source("code/functions_trajclass.R")
#source("code/subset_timeseries.R")
#noise_df <- get_data("data/simu/to_make_simu/noise_levels_dec.csv")
source("code/functions_trajclass.R")
source("code/subset_timeseries.R")
noise_df <- get_data("data/simu/to_make_simu/noise_levels_dec.csv")
#load("backup_dec.RData")

Abrupt trajectories

QuadraticA trajectories

QuadraticB trajectories

Linear trajectories

no changeA trajectories

no changeB trajectories

Calculate EWS metrics and DEV

Violin plots of DEV, for different categories and noise values

ploties <- list()
for (i in 1:length(noise_df$sr)) {
  for (largo in unique(EWS_data$length)) {
  
    p1<-ggbetweenstats(
      data = EWS_data[(EWS_data$error==noise_df$sr[i]&EWS_data$length==largo),],
      x = class,
      y = realEV,
      pairwise.display = "s", 
      centrality.plotting = F,
      xlab = "trajectory class",
      ylab = "DEV",
      title = paste0("DEV for different categories, length = ",largo,", demographic noise sr = ",
      noise_df$sr[i]),
    )
    ploties <- c(ploties,list(p1))
    
  }
}

Violin plots of AR1, for different categories and noise values

ploties <- list()
for (i in 1:length(noise_df$sr)) {
  for (largo in unique(EWS_data$length)) {
  
      p1<-ggbetweenstats(
      data = EWS_data[(EWS_data$error==noise_df$sr[i]&EWS_data$length==largo),],
      x = class,
      y = ar1,
      pairwise.display = "s", 
      centrality.plotting = F,
      xlab = "trajectory class",
      ylab = "Kendall tau for AR1",
      title = paste0("Trend in AR1 for different categories, length = ",largo,", demographic noise sr = ",
      noise_df$sr[i]),
    )
    ploties <- c(ploties,list(p1))
    
  }
}

Violin plots of SD, for different categories and noise values

ploties <- list()
for (i in 1:length(noise_df$sr)) {
    for (largo in unique(EWS_data$length)) {
  
      p1<-ggbetweenstats(
      data = EWS_data[(EWS_data$error==noise_df$sr[i]&EWS_data$length==largo),],
      x = class,
      y = sd,
      pairwise.display = "s", 
      centrality.plotting = F,
      xlab = "trajectory class",
      ylab = "Kendall tau for SD",
      title = paste0("Trend in SD for different categories, length = ",largo,", demographic noise sr = ",
      noise_df$sr[i]),
    )
    ploties <- c(ploties,list(p1))
    
    }
}

Violin plots of DEV (HALF-HALF method), for different categories and noise values

ploties <- list()
for (i in 1:length(noise_df$sr)) {
    for (largo in unique(EWS_data$length)) {
  
      p1<-ggbetweenstats(
      data = EWS_data[(EWS_data$error==noise_df$sr[i]&EWS_data$length==largo),],
      x = class,
      y = DEV.half,
      pairwise.display = "s", 
      centrality.plotting = F,
      xlab = "trajectory class",
      ylab = "DEV",
      title = paste0("DEV (HALF-HALF method) for different categories, length = ",largo,", demographic noise sr = ",
      noise_df$sr[i]),
    )
    ploties <- c(ploties,list(p1))
    
    }
}

Violin plots of AR1 (HALF-HALF method), for different categories and noise values

ploties <- list()
for (i in 1:length(noise_df$sr)) {
    for (largo in unique(EWS_data$length)) {
  
      p1<-ggbetweenstats(
      data = EWS_data[(EWS_data$error==noise_df$sr[i]&EWS_data$length==largo),],
      x = class,
      y = ar1.half,
      pairwise.display = "s", 
      centrality.plotting = F,
      xlab = "trajectory class",
      ylab = "log2 ratio AR1",
      title = paste0("Trend in AR1 (HALF-HALF method) for different categories, length = ",largo,", demographic noise sr = ",
      noise_df$sr[i]),
    )
    ploties <- c(ploties,list(p1))
    
    }
}

Violin plots of SD (HALF-HALF method), for different categories and noise values

ploties <- list()
for (i in 1:length(noise_df$sr)) {
    for (largo in unique(EWS_data$length)) {
  
      p1<-ggbetweenstats(
      data = EWS_data[(EWS_data$error==noise_df$sr[i]&EWS_data$length==largo),],
      x = class,
      y = sd.half,
      pairwise.display = "s", 
      centrality.plotting = F,
      xlab = "trajectory class",
      ylab = "log2 ratio SD",
      title = paste0("Trend in SD (HALF-HALF method) for different categories, length = ",largo,", demographic noise sr = ",
      noise_df$sr[i]),
    )
    ploties <- c(ploties,list(p1))
    
    }
}

Violin plots of AR1 (for the entire time series), for different categories and noise values

ploties <- list()
for (i in 1:length(noise_df$sr)) {
    for (largo in unique(EWS_data$length)) {
  
      p1<-ggbetweenstats(
      data = EWS_data[(EWS_data$error==noise_df$sr[i]&EWS_data$length==largo),],
      x = class,
      y = ar1.all,
      pairwise.display = "s", 
      centrality.plotting = F,
      xlab = "trajectory class",
      ylab = "AR1",
      title = paste0("AR1 for different categories, length = ",largo,", demographic noise sr = ",
      noise_df$sr[i]),
    )
    ploties <- c(ploties,list(p1))
    
    }
}

Violin plots of SD (for the entire time series), for different categories and noise values

ploties <- list()
for (i in 1:length(noise_df$sr)) {
    for (largo in unique(EWS_data$length)) {
  
      p1<-ggbetweenstats(
      data = EWS_data[(EWS_data$error==noise_df$sr[i]&EWS_data$length==largo),],
      x = class,
      y = sd.all,
      pairwise.display = "s", 
      centrality.plotting = F,
      xlab = "trajectory class",
      ylab = "SD",
      title = paste0("SD for different categories, length = ",largo,", demographic noise sr = ",
      noise_df$sr[i]),
    )
    ploties <- c(ploties,list(p1))
    
    }
}

EWSNet classification